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Abstract 

Quantum algorithms and complexity have recently been studied not only for discrete, but also 
for some numerical problems. Most attention has been paid so far to the integration problem, for 
which a speed-up is shown by quantum computers with respect to deterministic and randomized 
algorithms on a classical computer. In this paper we deal with the randomized and quantum 
complexity of initial- value problems. For this nonlinear problem, we show that both randomized 
and quantum algorithms yield a speed-up over deterministic algorithms. Upper bounds on the 
complexity in the randomized and quantum setting are shown by constructing algorithms with 
a suitable cost, where the construction is based on integral information. Lower bounds result 
from the respective bounds for the integration problem. 
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1 Introduction 



Potential advantages of quantum computing over deterministic or classical random- 
ized algorithms have been extensively studied by many authors for discrete problems, 



starting from Shor's paper on factorization of integers [14| and Grover's algorithm for 
searching databases Recently, a progress has also been achieved in quantum so- 
lution of numerical problems. The first paper dealing with the quantum complexity 
of a continuous problem was the work of Novak |13|, who established matching upper 
and lower bounds on the quantum complexity of integration of functions from Holder 
classes, based on the results on complexity of summation of real numbers from and 



1 11]. A general model of quantum computing for continuous problems has been devel- 
oped by Heinrich where the computation of a sum of real numbers is studied under 
various assumptions, and the results are applied to the integration problem. Another 
integration problem, computing path integrals, has been discussed in Recently, 
the approximation problem in discrete and continuous versions has also been treated in 
||5| and The linear problems of integration and approximation seem to be the only 
specific (and important) numerical problems discussed in the quantum setting so far. 
In the randomized setting, complexity results for problems such as integration, approx- 
imation or optimization are classical, see, e.g., ||l^ for an overview. 
In this paper we deal with the randomized and quantum solution of initial-value prob- 
lems. The complexity of this nonlinear problem was studied until now in the determin- 
istic worst-case and asymptotic settings, see, e.g., |§] and § for matching upper and 



lower bounds, or |10| for a discussion of complexity of initial-value problems on parallel 
computers. 

We show in this paper that a speed-up is achieved for initial-value problems by ran- 
domized and quantum algorithms over the deterministic ones. The results are summa- 
rized in Theorem 1, where we establish upper and lower bounds on the complexity of 
initial-value problems in the randomized and quantum settings. The upper bounds are 
obtained by defining algorithms based on the deterministic integral algorithm devel- 
oped in The procedure is shown allowing for an application of any algorithm for 
computing integrals (in the deterministic, randomized or quantum setting) to yield a 
new algorithm for initial-value problems in the respective setting. In the complexity 



analysis, the results on integration from [12| and []13| ] in the randomized and quantum 



settings are exploited. The comparison of the upper bounds on the randomized and 
quantum complexity to the worst-case complexity of initial- value problems shows that 
a speed-up is achieved in both non-deterministic settings. Lower bounds on the com- 
plexity refer to those for integration; the existing gap between the upper and lower 
bounds is discussed. 

2 Problem Formulation and Results 



We consider the solution of a system of ordinary differential equations with initial 
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conditions 

z'{t) = f{zit)), te[a,6], z{a)=rj, (1) 
where / : R'^ ^ R'^, z : [a, b] and ?? G R'^ (/(r/) / 0). 

Given an integer r > 0, /> G (0, 1], and positive numbers Dq, Di, . . . , Dr and H, we 
assume that the right-hand function / = [/^, . . . , /'^]^ belongs to the Holder class 

F-.P = {/:R'i^R'^| /gC7'-(R'^), [a7^(y)|<A, ^ = 0,l,...,r, 

<^l|y-^ir, y,^€R^ j = i,2,...,d}, (2) 

where 5*/-' represents all partial derivatives of order i of the jth component of /, and 
II • II denotes the maximum norm in R'^. To assure that / is a Lipschitz function, we 
assume that p = 1 for r = 0. 

We wish to compute a bounded function / on [a, b] that approximates the solution z. 
Letting {xi} be the uniform partition of [a,b], Xi = a + ih with h = (b — a)/n, the 
function I will be produced by an algorithm cj), based on approximations to z{xi), 
i = 0, 1, . . . , re. 

We now discuss the error and complexity models in the worst-case deterministic, ran- 
domized and quantum settings. In the worst-case deterministic setting, the error of (p 
at / for the problem (|^) is defined by 

e{^,f) = sup ||z(t)-/(t)||, 
and the error in the class F^'f by 

gWorst(^^^r,p)^ sup e(,/.,/). (3) 

We assume that the values of / or its partial derivatives can be computed at given points 
by a subroutine. The cost of an algorithm (j) is measured by a number of subroutine calls. 
For a given e > 0, by the e-complexity of the problem, comp™°'''^*(i^''''', e), we mean the 
minimal number of subroutine calls (taken among all possible algorithms) sufficient to 
solve the problem with error at most e, i.e., the minimal cost of an algorithm (p taken 
among all </> such that e"™"*((/), F^'^p) < e. 

In the randomized setting, we allow a random selection of points at which the function / 
is evaluated, so that the output of an algorithm is a random variable (on a probability 
space (ri, S, P)). Let the mappings uj ^ ^ of (/) be random variables for each 
/ E F'^'P By an algorithm in the randomized setting, we mean a tuple 

<^ = ({ao(-),<(-),---,«n(-)Wf^,^), (4) 
where ip \s a, mapping that produces a bounded function 

r(t)=V«(/),a^(/),...,a-(/))(t), (5) 
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t S [a, b], based on ag (/), (/)> • • • , o-nif)- The error of at / for the problem (||) is 
defined by 

e-(0,/)= sup \\z{t)-l'^{t)\\. (6) 

te[a,b] 

We assume that the mapping G — > e'^(i?!>, /) is a random variable with values in R, 
for each / G F'^'^. The error of (p in the class F'^'^ is given by 

^rand^^^^r.p)^ SUp (Ee^ (</., V2 ^ (7) 

where E is the expectation. 

As in the worst-case setting, we measure the cost of an algorithm by a number of 
subroutine calls that are needed to compute an approximation. For a given e > 0, by 
the e-complexity of the problem, comp'''^'^'^(F''''', e), we mean the minimal cost of an 
algorithm taken among all (j) such that e''^'^'^(</>, F'''^) < e. 

In the quantum setting, the output of an algorithm is a random variable (taking a 
finite number of values), but the reason of randomness is different than that in the 
randomized setting. On a quantum computer, where, roughly speaking, basic objects 
are qubits (elements of a two-dimensional complex space Hi) and allowed operations 
are unitary transformations of the the tensor product of a number of copies of Hi, 
randomness is a result of quantum measurement operations. For a detailed description 
of the framework of numerical quantum computing one is referred to Q, where the 
notions of quantum measurement, quantum query, quantum algorithm and complexity 
are defined and thoroughly discussed, and applications to summation and integration 
problems are studied. For a condensed discussion of randomized and quantum settings, 
in particular for the integration problem, one is also referred to 0. 
By a quantum algorithm for solving our problem we mean a tuple (Q), where af{f) 
are random approximations, in the quantum sense, to z{xi) for each /. The error of 4> 
at / is defined by (P). 

Let < 6 < 1/2. The error of (p in F^'P in the quantum setting is defined @ by 

gquant^^^ pr,p ^ ^ -^^f | ^[ p| e-(,^, /) > a } <6 }. (8) 

Note that for a given e > the bound e'^{(j), /) < e holds with probability at least 1 — 6 
for each / iff e<i"^'^*((/>, F'^'P, 5) < e. 

In the quantum setting the value of 5 is usually set to 5 = 1/4. The error probability 
can then be reduced to any 5 by computing (componentwise) the median of c log 1 /5 
repetitions of the algorithm, where c is a positive number independent of 5, see 
Lemma 3. For our problem, the procedure of increasing the probability of success 
can be applied at different levels which influences a logarithmic part of the cost of an 
algorithm, so that we shall describe it in more detail and discuss after the proof of 
Theorem 1. 

On a quantum computer, the right-hand side function / can be accessed through a 
query that returns, for a given point, a value of a component of /. Roughly speaking. 
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a query on a class of real functions is defined as a transformation Q that associates 
with each function p a unitary mapping Qp defined on a Hilbert quantum space. For 
a detailed discussion of what is meant by "returning a value" of a function, and how 
a query is implemented in the quantum setting, the reader is referred to Q or 
The cost of an algorithm (p is measured by a number of quantum queries that are 
needed to compute an approximation. (In upper bounds in Theorem 1, classical evalu- 
tions of / or its partial derivatives are also taken into account.) For a given e > 0, by 
the quantum e-complexity of the problem, comp'i^^°*(i^'"''', e, 5), we mean the minimal 
cost of a quantum algorithm (p taken among all (j) such that e"'^^^^{(j), F'^'P, 6) < e . 
We prove in this paper upper and lower bounds on the randomized and quantum com- 
plexity of initial- value problems (|l]). The upper bounds, summarized in the following 
theorem, will be next compared to the known lower bounds on deterministic complexity 
to show that a speed-up is achieved in both settings. Lower bounds in the randomized 
and quantum settings are derived from a simple argument in the case d > 2, and are 
also included. We take below log = log2 (although the base of the logarithm is not 
crucial). 

Theorem 1 For the problem we have that 

( r+p+3/2 \ 
f 1\ (r+p+l/2)(r + p+l) 1\ 

/ r+p+2 \ 

comp^-'^*(F'-''', 8,6) =0 I (log i + log i) . (10) 



For d>2, 



and, forO<6< 1/4, 



compi^^^*(F'''^e,(5) > comp^"'^'^*(F^'^ e, 1/4) = n ^^Qy^'^'^ . (12) 

The constants in the "O" and "Q" notation only depend on the class F^'P , and are 
independent of e and 5. 

Upper bounds (P) and ( p!o|) will be derived by defining suitable algorithms, while the 
lower bounds (0) and (^) are equal to those on the complexity of randomized or 
quantum computation of integrals of a function of one variable. 

Before giving the proof, we make some comments on these results. If the values of / 
or its partial derivatives can only be accessed, the deterministic worst-case complexity 
of the problem (|l|) is of the order e~i/(^+/')j see Theorem 3 in the next section. Since 

^ + ^ + ' '-±J^ ,<^, (13) 



{r + p + iy {r + p+l/2){r + p+l) r + p 
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both randomized and quantum computation yield an improvement over the determin- 
istic setting over the entire range of r and p (we neglect the logarithmic factors). For 
instance, if r = and p = 1, the worst-case complexity in the deterministic setting is 
of the order e~^, in the randomized setting it is bounded from above by e~^/^, while 
on the quantum computer by 



The lower bounds coincide with those for the integration problem, see (16) and (|T] 
To see what the size of the gap between the bounds is, note that the reciprocal of the 
exponent in 1/e in the quantum case is such that 

lim f (^ + ^ + ^)' -(r + p)Uo, 

so that it behaves for large r as r -|- p, while the lower bound depends on r + p + 1. In 
the randomized setting, the reciprocal of the exponent in 1/e behaves for large r like 
r + p, while the lower bound depends on r + p + 1/2. 

We now recall results on randomized and quantum computation of integrals, as well as 
those on deterministic solution of initial-value problems. 

3 Randomized and Quantum Computation of Integrals 
and Deterministic Solution of Initial- Value Problems 



Quantum complexity of integration has been first studied by Novak [|1^]. The problem 
is to approximate the integral 



I{g)= J g{x)dx (14) 

[0,1]« 

for functions g : [0, 1]* — > R from a Holder class with r>0, 0<p<l 

F'^f' = {geC'\[0,ir) I \9iy)\<Do, giy) - giz)\ < H \\y - z\\p , y,zG [0,1]^}, 

(15) 

where d'^g represents all partial derivatives of order r of g. 

The error of an algorithm cj) at g for the integration problem (14) in the randomized 
and quantum settings is defined by 

e'^{cP,g) = \I{g)-A-{g)l 

where A'^(g) is the output of (j) (in the worst-case setting the definition is the same, only 
the output is deterministic). The other definitions of errors in the class of functions 
and complexity remain the same as for the problem (|l]), with F^'^ replaced by F^'^. 
Based on the results on the computation of the mean of n numbers given by Brassard 



et al ||l| (upper bound), and Nayak and Wu |Tl|] (lower bound), Novak [13| showed 



the following result in the quantum setting. For the result in the randomized setting. 



see IjlJ, p. 62. Let 7 = (r + p)/s. 
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Theorem 2 ( [jl^ , For the problem ^14) we have 

comp'^^°^(F"''', e) x £-1/(7+1/2) ^ (^g) 

^.Q^pquant(^r,p^g^ 1/4) X £"1/(^+1). (17) 



An upper bound in (|16| ) can be achieved by random algorithms with a finite number 
of output values. 

Consider now the solution of initial- value problem (Q) in the deterministic setting. This 
problem has been considered in a number of papers, see, e.g., ||8| or [^J. The following 
result is a straightforward modification of Corollary 4.1 from The modification is 
needed, since the class of functions considered in § , consisting of r times continuously 
differentiable functions with bounded derivatives, is to be replaced with the Holder class 
F'^'^. For the modification in lower bounds, one is referred to the proof of Theorem 
3.1 from Ip, where the functions in the construction must be replaced by suitable 
functions from the class F'^'f'. The upper bounds will be derived again in the sequel, as 
a by-product in the proof of Theorem 1. 

Theorem 3 ([^) In the deterministic setting, the complexity of (0J satisfies: 

- if the values of f or its partial derivatives are only accessible, then 

comp'"°''\F'^'P, e) - £-i/(^+'') , (18) 

- if arbitrary linear functionals are accessible, then 

com.-p"'"'\F'''f, e) x £-i/(r+P+i) _ (^g) 

■ 

(The lower bound in ( [T9| ) holds true not only for linear functionals, but also for a class 
of nonlinear functionals, see (P.) Relation ( p!9|) will play an important role in the proof 



of Theorem 1, while (18) will serve as a point of reference in evaluating a speed-up 
obtained due to randomization or due to quantum computations. 
Let us now recall the algorithm that leads to the upper bound in (|l|). It requires the 
computation of integrals of /, and is defined as follows. 

Take = rj. Given y* {y* = z{xi)), we let zl{t) be the solution of the problem 

z'{t) = f{z{t)), te[x^,Xi+l], z{xi)=y*, (20) 

and we set 



r+l -j^ 



= E it - ^i)' . t e \xi,Xi+i] . (21) 



j=0 



Then we define 



yhi = y:+ / fimt))dt, (22) 
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i = 0, 1, . . . , n — 1, and finally 

l{t)=l*{t) ^OVte[x^,X^+l]. (23) 

(The function / is piecewise continuous. It is also possible to define it to be continuous 
on [a, b].) 

4 Randomized and Quantum Solution of Initial- Value 
Problems 

We now define randomized and quantum algorithms for the solution of (||). Let w* be 
a polynomial 

<(y) = E V^'^(2/:)(y-y*y, (24) 

where f^^\y*)z^ is meant to be the value of the j-linear operator f^^\y*) at {z,z . . . ,z) 
(j times). The values of w* can be computed through evaluation of partial derivatives 
of components of / of order 0,1, ... ,r. Equality (22) can be equivalently written as 

^i+i 1 

yhi = y:+ J w*{r*{t)) dt + /i^+^+i J g,{u) du , (25) 

Xi 

where 

for u € [0, 1] and z = 0, l,...,n — 1 . One can verify that gi belongs to C^'''\[0, 1]), the 
derivatives of gi are bounded by constants that depend only on the parameters of the 
class F'^'P (and are independent of i, y* and h), and 

\\gl'\u)-g^/-\u)\\<H\u-u\P, 

u,u £ [0, 1], for some constant H depending on the parameters as above. 

The algorithm (deterministic, randomized or quantum) for solving (||) is defined as 

follows (we omit the argument uj in random variables). Let ao(/) = yo = V- Given 

c-iif) = yi) we consider functions gi defined by ( ^6|) for yi (that is, the polynomials 

I* and w^ based on y* are replaced by the polynomials k and Wi based on y^), and 

1 

compute some approximations Ai{f) to the integrals / gi{u) du in (|25|). The algorithm 



is defined by setting aj+i(/) = yi+i, where 

yi+i = yi+ I Wi{k{t))dt + h^+P+^Ai{f), (27) 



and the approximation on is given by 

l{t) = hit) , (28) 
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i = 0, 1, . . . , n — 1. Finally, we set ip{aQ{f), . . . , an{f)){t) = l{t) for t G [a, b]. 
The approximations Ai{f) may be obtained by deterministic, randomized or quantum 
algorithms. In the randomized and quantum settings, we demand random variables 
Ai{f) to satisfy 



gi{u) du-Ai{f) 



< El with probability at least (1 — 6) 



l/n 



(29) 



for i = 0, 1,... ,n — 1 (and all in), for some ei. It will be shown later on that a 
satisfactory choice is ei = h. 

For illustration, we specify the algorithm above in the case r = 0. It may be considered 
as a modification of Euler's method and is defined as follows. Given yi, we compute 

yi+i=yi + hf{yi) + h^+PAi{f), 



where Ai{f) is a (deterministic, randomized or quantum) approximation to {l/h^) J{f{yi+ 



uhf{yi)) — f{yi)) du with error at most ei, and probability at least (1 — 5)^/^ (in non- 
deterministic cases), i = 0,l,...,n — 1. The approximation to z = z{t) is defined on 
[xi,Xi+i] by l{t) =yi + f{yi){t - Xi). 

Proof of Theorem 1 

Consider first the quantum setting. We use the quantum algorithm of Novak to compute 
Ai{f) (componentwise in a statistically independent way). Due to (^) with s = 1, we 
have that for each i the inequality 



(30) 



gi{u)du - Ai{f) 

holds with probability at least 3/4, and the quantum query cost 0{e-^^'^''^'^^^). (Since 
Qi has d components, only the constant in the "O" notation is different than that in (p!7|), 
which is a result of computing the median of a suitable number of repetitions to increase 
the probability of success in each component to (3/4)^/"^.) Taking componentwise the 
median of k repetitions, 



= e log 



1 



1-(1-(5)V"^ 

we arrive at an approximation Ai{f) (the same symbol is used to denote this new 



approximation) such that ( ^91) is satisfied. Since 



log 



1 - (1 - 5)1^ 



< c(log n + log 1/6) 
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(where c is independent of n and 5), the cost of computing Ai{f) is of order 0{ei (log n+ 

logl/(5). Thus, with the cost 0(5]^ ^^^''^''^^^ n (log n + log 1/(5)) we assure that the 
bounds 



gi{u) du - Ai{f) 



< El for z = 0, 1, . . . , n — 1 



(31) 



hold simultanously with probability at least 1 — 6. The total cost of the algorithm 
additionally includes 0{n) classical subroutine calls necessary in the deterministic part 
of the algorithm. 

Let Cj = z{xi) — iji. Since the solution of (|l|) satisfies 



z{xi+i) = z{xi) + J f{z{t))dt, 

we get from (^) that 

ei+i = ei +J\f{z{t))-f{kmdt 

+ Yifihit)) - wiikm dt - h'+p+^Mf) . 



(32) 



(33) 



We shall derive from ( |33D a difference inequality for 1 1 ej 1 1 . Let Zi denote the solution of 
the local problem (pH) with the initial condition Zi{xi) = yi. By triangle inequality, 



\\f{z{t)) - fum\ < wfizit)) - fim)\\ + wfim) - fmm- 



(34) 



To estimate the first term, we note that the dependence of the solution on initial 
condition yields that 



\z{t) - Zi{t)\\ < exp{Lh)\\z{xi) - yi\\, t £ [xi,Xi+i], 



(35) 



where L is the Lipschitz constant of /. Writing the remainder of Taylor's formula in 
the form 



m - hit) = J {4"^'\et + (1 - e)x,) - +'^(xi))(t - xiY+H^ - oy/n de , 



t G [xj, Xj+i], and checking that z\''^^^ is a Holder function with exponent p (and a 
constant that only depends on the parameters of the class F''''P and is independent of 
i,yi,n), we arrive at 

\\zi{t) - k{t)\\ < M h'+P+\ (36) 

j+i], where the constant M only depends on the parameters of the class F^'P. 
The last two inequalities together with (|33D, (34) and (^) yield the relation (satisfied 
with probability at least 1 — 6) 



\ei+i\\ < \\ei\\ (1 + hLexp{Lh)) + LMW+p+^ + N'+p+^Ei 



(37) 
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0,l,...,n — 1. Take now £i = h. By solving difference inequality (p7|) , we get that 

max lleill =0(/i^+''+i) (38) 

0<i<n 



with probability at least 1 — 6. Finally, for any t G we have due to (pq) , (|3q) 

and (ID that 

\\z{t) - i{t)\\ = Mt) - km < \\z{t) - zi{t)\\ + \\m - km = oih-+p+'), m 

with the constant depending only on the parameters of the class F'^'P, and probability 
at least 1 — 6. 

The quantum query cost of the considered algorithm is of order 0{n^'^~^P~^'^^^ (r+p+i) Q^^g 
log 1/6)). For e > 0, we now take the minimal n such that the upper bound in (|3^ ) 
does not exceed e, n x (l/e)i/('-+p+i) , and the desired upper bound on the complexity 
(0) follows. 

In the randomized setting we use ([T6|) . Take an algorithm (j) for approximating integrals 
with s = 1 such that e^^^'^{(f), F^'P) < ei/2. By the Chebyshev inequality we have that 



1 

g{x)dx-A^{g) 







> £1 ^ < 1/4, 



for any g € F^'f, where A'^(g) is the output of the algorithm. As in the quantum 
case, the inequality ( |30[ ) holds with probability at least 3/4, but the cost is now 
0{£i ^/^'^'^f'^^/'^^)^ Proceeding further on similarly as in the quantum case, one gets that 
the bound (H) holds with probability at least 1—6, and cost 0(n('"+''+^/^)/(^+''+^/^) (log n+ 
log 1/6)). We denote the resulting random algorithm for initial value problems by (p. 
Take now n to be the minimal number for which the upper bound in (p9|) is at most 
e/2, n X (l/e)i/(''+/'+i) , and observe that 

Ee-(0,/)2= J e^{^J)^dP{u;)+ J e'^ {^J)^ dP{u;) < K^6 + /A 

e-(^,/)>e/2 e-(^,/)<£/2 

for all / € F'^'P, where X is a positive constant that depends only on the parameters 
of the class F^'f such that e'^{(j), f)<K (such a constant exists, since we may assume 
that I ^'^((7) I < 2Dq; otherwise A'^^g) = would be a better approximation). Hence, 
the algorithm cj) with 6 = 3e^/4K^ satisfies e'''^°'^((^, F^'f) < e. Looking at the cost of cj) 
we see that the upper bound (P) is proven. 

Consider now the deterministic setting. The known upper bounds on the complexity 

can be derived again as follows. If the values of / or its partial derivatives can only 

1 

be accessed, then we are able to approximate / gi{u) du within the error £1 with cost 



0(e^ i/(''+p)^_ Since n integrals are to be approximated, the total cost of the algorithm 
(for £1 = h), with the error bound (PD, is of order 0{n ■ v}/'^''+p^). This proves the 
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upper bound in (p!8|). If exact computation of the integrals is allowed in the model, 
then we simply take ei = 0, which leads to the upper bound in (19). 




We now pass to lower bounds. Consider the quantum setting. Let g G F^''^ with r > 1 
and /9 G (0, 1], or r = and p = 1, and with s = 1. Consider a two-dimensional problem 

= giuit)), tG[0,l], 

with initial conditions n(0) = 0, v{0) = 0. The function g can be extended to R such 

that the right-hand side function in ( ^0| ) belongs to the class F^'f with suitably chosen 

t 

parameters. The solution of (|^) is given by u{t) = t and v{t) = J g{s) ds. 



Let e > and let (j) be any quantum algorithm for solving ([l|) with error e°'^^'^^{(j)^ F'^'P , 1/4) < 

e, and quantum query cost c{e). When applied to ([40|), the algorithm (f) gives an ap- 

1 

proximation to f (1) = / g{t) dt with error at most e, with probability at least 3/4. Due 



to the lower bound for integration in (|l^) the cost must be at least of order 
queries on g, which yields that c{e) = Q.{£~^/^''^'^f^'^'>), and proves lower bound (12). 



In the randomized setting we use similar arguments to show (|11|), adjusted to the error 
formula (0). Since the same arguments apply in both settings for any d > 2, the proof 
of Theorem 1 is completed. 

Let us finally note that the logarithmic factor in the upper bound (|lO|) depends on what 
stage the median is computed at in the algorithm. If we ask the inequality (^) to hold 
with probability at least (3/4)"^/", proceed up to the final step (^) with 5 = 1/4, and 
after that compute the median of clog 1/5 repetitions of the entire algorithm, then the 
factor log 1/e + log 1/5 in (|lO|) would be replaced by log 1/e • log 1/5. 
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